Mechanism for Spontaneous Growth of Nanopillar Arrays 
in Ultrathin Films Subject to a Thermal Gradient 
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Several groups have reported spontaneous formation of periodic pillar-like arrays in molten poly- 
mer nanofilms confined within closely spaced substrates maintained at different temperatures. These 
formations have been attributed to a radiation pressure instability caused by acoustic phonons. In 
this work, we demonstrate how variations in the thermocapillary stress along the nanofilm interface 
can produce significant periodic protrusions in any viscous film no matter how small the initial trans- 
verse thermal gradient. The linear stability analysis of the interface evolution equation explores an 
extreme limit of Benard-Marangoni flow peculiar to films of nanoscale dimensions in which hydro- 
static forces are altogether absent and deformation amplitudes are small in comparison to the pillar 
spacing. Finite element simulations of the full nonlinear equation are also used to examine the array 
pitch and growth rates beyond the linear regime. Inspection of the Lyapunov free energy as a func- 
tion of time confirms that in contrast to typical cellular instabilities in macroscopically thick films, 
pillar-like elongations are energetically preferred in nanofilms. Provided there occurs no dewetting 
during film deformation, it is shown that fluid elongations continue to grow until contact with the 
cooler substrate is achieved. Identification of the mechanism responsible for this phenomenon may 
facilitate fabrication of extended arrays for nanoscale optical, photonic and biological applications. 



I. INTRODUCTION 



The manufacture of ultra small optical and electronic 
components is nowadays based on optical lithography 
techniques whereby a geometric pattern defined by a pho- 
tomask is transferred onto a photosensitive resist layer 
by exposure to UV light. Various chemical treatments 
are then used to embed the positive or negative image 
of this pattern onto a material film beneath the pho- 
toresist. While this commercial technique can generate 
feature sizes below 100 nm, there are certain disadvan- 
tages inherent in the patterning process [1]. For example, 
multiple step-and-repeat processes are required for de- 
position, exposure and removal of the photoresist layers 
for constructing three dimensional components. Inhomo- 
geneities in the photoresist layer thickness, composition, 
exposure dose or developer concentration can cause sig- 
nificant surface roughness and scattering losses which di- 
minish performance of optical or electronic components. 
Optical lithography is also inherently a two-dimensional 
technique whereby three dimensional components are 
fabricated layer upon layer. The process requires that 
the supporting substrates be rigid and fiat, posing chal- 
lenges for the fabrication of curved or complex shaped 
components. In an effort to eliminate such constraints 
while reducing fabrication time and cost, researchers have 
been exploring alternative, lower resolution patterning 
techniques such as ink-jetting [5], gravure printing [3], 
direct-write [1], micro- moulding [S] and nanoimprinting 
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[BHS] . These methods are more adaptable to new materi- 
als and pattern layouts; however, multiple etching steps 
are still required and device performance is still not com- 
parable to those fabricated by conventional means. The 
materials of choice tend to be inks, colloidal suspensions 
and polymer melts [5], which are not only less costly but 
whose composition can be tuned to optimize functional- 
ity. 

Some groups have been investigating less conventional 
means of film patterning by exploiting the self-assembling 
character of structures formed by hydrodynamic insta- 
bilities in thin films. Examples include dewetting in- 
duced by chemically templated substrates [TD], capillary 
breakup on rippled substrates [TT], island formation in 
ferroelectric oxide films [12], elastic contact instabilities 
in hydrogels |13| and evaporative instabilities in metal 
precursor suspensions [2] . The use of fluid instabilities 
for controlled formation of large area, periodic arrays pro- 
vides an interesting approach for future development of 
non-contact, resistless lithography. 

It is well known that liquid films with dimensions in the 
micron to nanometer range manifest exceedingly large 
surface to volume ratios. As such, small liquid structures 
can respond instantaneously to external modulation of 
surface forces. This sensitivity to surface manipulation 
has been successfully used to control the motion of small 
liquid volumes for micro-, bio- and optofluidic applica- 
tions [H]. For example, tangential stresses based on ther- 
mocapillary forces have been used to steer [TBHT5] . mix 
[TOl and shape "M. thin films and droplets on demand. 
Since the surface tension of liquids varies with tempera- 
ture, thermal distributions can be applied directly to a 
supporting substrate to generate lateral gradients which 



drive the flow of liquid toward selected regions of a sub- 
strate. In this work, we examine systems comprising of 
liquid nanofilms subject to a transverse temperature gra- 
dient which have been observed to produce nanopillar ar- 
rays which grow and elongate in the direction of a cooler 
target substrate. The spontaneous formation of 3D large 
area arrays offers exciting possibilities for non-contact, 
resistless, one step fabrication of optical and photonic 
structures. Since solidification of the emergent molten 
structures occurs in-situ upon removal of the thermal 
gradient, it is anticipated that the resulting nanostruc- 
tures will manifest specularly smooth interfaces, a dis- 
tinct advantage for optical applications. 



Cooler substrate r> 




A. Formation of nanopillar arrays in molten 
polymer nanofilms 

The typical experimental setup leading to spontaneous 
formation of nanopillar arrays is shown in Fig. llTa) . Poly- 
mers such as polystyrene (PS) or poly(metylmetacrylate) 
(PMMA) are first spun cast onto a clean, flat silicon wafer 
to an initial thickness ho of the order of a few hundred 
nanometers. The coated wafer is then overlay with a 
second silicon wafer containing vertical spacers along the 
periphery to ensure an air gap above the polymer film. 
The wafer separation distance, do-, is normally several 
hundred nanometers. The bottom and top wafers are 
maintained at different temperatures above the polymer 
glass transition temperature to ensure a fiowing liquid 
film. In all the experiments reported in the literature, 
AT = T2 - Ti « 10 - 50° C. Next, we review the ex- 
perimental results of three independent groups reporting 
observations and measurements of nanopillars arrays. 



1. Experiments by Chou et al. 

Chou et al. [HI [21] appear to have been the first group 
to report nanopillar formation in ultrathin polymer films. 
In their experiments, they studied low molecular weight 
PMMA (approx 2K), which was first spun cast to a film 
thickness of 100 nm onto a cleaned silicon wafer and then 
annealed at 80 °C to drive off residual solvent. The an- 
nealed film was then placed within the assembly shown 
in Fig. [l[a) , where the top wafer had been treated with a 
nonstick coating to prevent polymer attachment after so- 
lidification. The underside of the top wafer was either flat 
or patterned with a rectangular relief structure a few tens 
of microns in width and about 0.3/im tall. In all exper- 
iments reported, there was no imposed temperature dif- 
ference between the top and bottom wafers {T2 — T1 = 0). 
Instead, the entire assembly was cyclically heated from 
room temperature to either 130 °C or 170 °C, well above 
the polymer glass transition temperature Tg — 103 °C 
[25] to ensure a softened film. The heating cycle per- 
sisted for 5-80 minutes with no noticeable difference in 
pattern formation if the air gap was replaced by a vac- 



FIG. 1: (a) Sketch of experimental setup for formation of 
nanopillar arrays. Initial thickness of flat nanofilm is denoted 
by ho; gap spacing in between silicon substrates is denoted 
by do- Length scale Amax represents theoretical prediction 
for pillar spacing; \°^p^ represents experimentally measured 
values, (b) AFM image of PMMA pillars [21]: ho = 95 nm, 
do = 260 nm, AT unknown, A'^''^* = 3.4 /xm. (c) Optical 
micrograph of PS pillars [2^: ho — 100 nm, do = 285 nm, 
AT = 46 °C, A™P* = 2.9 ± 0.6 ^m. (d) AFM image of PMMA 
piUars 23/. ho = 100 nm, do = 163 nm, AT = 10 °C, A™^' = 
6.5 ^m. 



uum at 0.3 Torr. In cases where the PMMA coated wafer 
was not overlay by a top wafer and simply exposed to 
open air, no protrusions were observed to form. When 
the top wafer was placed in close proximity to the melt 
surface i.e. {do — ho) ~ 165 nm , nanopillar arrays with 
in-plane hexagonal symmetry were obtained, as in the 
image shown in Fig. [l[b). These elongations were mea- 
sured to have a diameter and pitch (i.e. pillary spacing) 
of a few microns; their overall height closely matched the 
gap distance do separating the two wafers. AFM images 
of the resulting structures after solification revealed pil- 
lars with a flat top and fairly straight sidewalls. Chou 
et al. attributed the formation of these elongations to 
an image-charge induced electrohydrodynamic instabil- 
ity caused by non-uniform distribution of charges on the 
relief surface. Chou et al. also noted that thermal gradi- 
ents might be playing a role but that Rayleigh-Benard or 
Benard-Marangoni cellular convection was unlikely since 
the initial film thicknesses were far too small to overcome 
the relevant critical numbers required for instability [H] . 



2. Experiments and m,odeling efforts by Schaffer et al. 

Soon thereafter, Schaffer and co-workers [52 [SS] [27] 
used a similar setup as in Fig.ma) where the two confin- 
ing wafers were purposely set to different temperatures 
such that T2 > Ti. They first spun cast high molecu- 
lar weight films of PS (Tg = 95 °C [Sg, mol. wt. 108 
kg/mol) dissolved in toluene onto a silicon wafer down 
to an initial thickness 80 nm ^ ho ^ 130 nm. It appears 
that these films were not annealed to drive out residual 



solvent after spin casting, which may have led to overes- 
timates in the reported values of ho (discussed further in 
Section III). The wafer separation distance ranged from 
100 nm < do ^ 600 nm. The bottom wafer was then 
heated to T2 — 170 °C; the top wafer was cooled to a 
temperature above Tg such that AT = T2 — Ti ranged 
from 10 < AT < 55 °C. The small wafer separation 
distances give rise to very large transverse thermal gra- 
dients of the order of AT/do ~ 10^ - 10^ °C/cm. After 
subjecting the PS film to the thermal gradient overnight, 
the sample was quenched to room temperature and the 
top wafer removed. As in Chou et al. , the top wafer had 
been treated with a silanized monolayer in order to pre- 
vent adhesion of the PS. After solidification and removal 
of the top wafer, the films were observed to contain pe- 
riodic nanopillar arrays, as shown in Fig.fHc). To study 
the influence of the wafer separation distance do on the 
pillar formation process, Schaffer et al. used a tilted plate 
geometry in all their experiments where the top wafer 
was inclined with respect to the bottom one by about 
1/im over a distance of 1 cm, corresponding to an incli- 
nation angle of about 0.0057°. This modification allowed 
simultaneous measurement of the array pitch as a func- 
tion of do within a single run. Schaffer et al. conducted 
a comprehensive set of experiments and determined the 
influence of the initial film thickness ho, the wafer separa- 
tion distance do, and temperature drop AT on the pillar 
separation distance A. They ruled out any electrostatic 
effects by purposely grounding the confining wafers. 

As noted both by Chou et al. [2T] and Schaffer et 
al. [26], films ranging in thickness from millimeters to 
centimeters subject to a transverse thermal gradient are 
known to develop cellular instabilities which lead to pe- 
riodic surface deflections at the air/liquid interface due 
either to Rayleigh-Benard (RB) or Benard-Marangoni 
(BM) convection [23. These instabilities, however, gen- 
erate very shallow corrugations and not pillar-like pro- 
trusions as observed in nanofilms. Onset of instability 
requires that the critical Rayleigh number -Rflonsot for 
buoyancy driven flow (which scales as /i^) or the criti- 
cal Marangoni number Maonact for thermocapillary flow 
(which scales as h"^) exceed 660 - 1700 or 50-80, respec- 
tively, depending on the boundary conditions. In the 
nanofllm experiments, the corresponding values are es- 
timated to be Ra w 10~^^ and Ma w 10^^, orders of 
magnitude less than required for onset of instability. 

Schaffer et al. therefore proposed a different mecha- 
nism for instability and interfacial deformation based on 
a novel radiation pressure model. They hypothesized 
that low frequency acoustic phonons (AP) can reflect 
coherently from the interfaces of the molten film over 
distances of the order of the film thickness despite that 
the melt is in an amorphous state. These low frequency 
modes are postulated to generate a significant desta- 
bilizing radiation pressure while conducting little heat. 
By contrast, the high frequency modes are expected to 
propagate diffusively with little interfacial resistance and 
therefore little interfacial pressure. These modes, how- 



ever, are essential for establishing the steady-state heat 
fiux across the air and melt layers. The mechanism de- 
scribed represents a kind of acoustic analogue of the ra- 
diation pressure caused by optical phonon reflections in 
closely spaced metal plates placed in vacuum, known to 
generate the Casimir interaction force [29]. Since the 
air/melt interface is liquid-like and therefore deformable, 
the acoustic phonons in the polymer melt are believed 
to generate an outwardly oriented radiation pressure, 
which counteracts the stabilizing force of surface tension; 
infinitesimal surface defiections can therefore grow into 
sizeable protrusions. Schaffer et al. developed a detailed 
hydrodynamic model based on the slender gap approxi- 
mation for describing the evolution equation for the film 
thickness, h{x,y,t). A linear stability analysis of this 
evolution equation leads to an analytic expression for the 
wavelength corresponding to the fastest growing unstable 
mode, namely 
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where 7 denotes the surface tension of the polymer melt. 
Up is the speed of sound in the polymer melt and n = 
kair/kmeit dcnotes the ratio of thermal conductivity of 
air to that of the polymer melt. The superscript AP dif- 
ferentiates this expression from the one to be derived for 
a thermocapillary model (TC). The material constants in 
Eq. (IT]) are evaluated at the substrate temperature T2. 
The parameter Q represents the acoustic quality factor 
determined from the phonon reflection and transmission 
coefficients corresponding to the four media constituting 
the system, namely the bottom silicon wafer, the poly- 
mer melt, the overlying air layer and the top silicon wafer. 
Positive values of Q lead to fllm destabilization and the 
formation of nanopillar arrays. Schaffer et al. compared 
the prediction for A^^x directly with the pillar spacings 
obtained in experiment, A'^''^*. A least squared fit of the 
experimental data to the model with Q and Up as fitting 
parameters produced good agreement (see dashed curves 
in FigMb) ). In particular, it was shown that the value 
of Q did not vary with ho, do or AT. The acoustic qual- 
ity factor Q seemed to depend on the choice of substrate; 
Q = 6.2 was obtained for the silicon/air/PS/silicon sys- 
tem, while Q = 83 for films supported by a silicon wafer 
coated with a 100 nm layer of gold. Unfortunately, it was 
reported [17] that the measurements of X'^^p^ included 
not only pillar formations but lamellar structures, spirals 
and other periodic formations caused either by defects in 
the initial film or by prolonged contact with the cooler 
substrate. In many cases, protrusions had undergone re- 
organization while in contact with the cooler wafer. In 
addition, measurements of pattern periodicity were ob- 
tained long after contact with the upper substrate and 
subsequent solidification. Comparison of these measure- 
ments to a model based on linear instability is therefore 
problematic. 

Schaffer et al. concluded that they had uncovered a 
novel instability in nanofilms induced by a radiation pres- 



sure from interfacial reflections of low frequency acoustic 
phonons. They noted that the frequency dependence for 
propagation of acoustic phonons with large mean free 
path is highly unlikely in low molecular weight polymers 
and that the instability would not be observed in such 
systems since they lack the necessary glassy rheological 
response [35]. In a separate study, Schaffer et al. [50] 
also conducted experiments with relief structures pat- 
terned with complex patterns held in close proximity to 
the polymer melt interface. The smallest values of do /ho 
lead to well defined replicas in the polymer film. 



3. Experiments by Peng et al. 

Shortly following the work of Schaffer et al. , Peng and 
co-workers [23] used a similar assembly as in Fig. [l]^a) 
to study PMMA films with ho ~ 100 nm, T2 = 160 °C, 
130°C < Ti < 150°C and 110 < do < 210nm. They were 
able to obtain nanopillar arrays after about 0.5 — 2.5 hrs; 
however, they did not conduct a parametric study nor 
compare their measurements of the pillar spacing with 
the prediction of Schaffer et al. . Fourier transforms of 
the nanopillar arrays showed well defined hexagonal sym- 
metry in some cases, as shown in Fig. flld). In other ex- 
periments, the pillar formations adopted either stripe or 
spiral symmetry. Peng et al. used a simple energy min- 
imization argument first introduced by Schaffer et al. to 
show that pattern selection between stripe and hexago- 
nal arrangements is merely controlled by the thickness of 
the overlying air film, while spiral formations are likely 
caused by point defects in the film. In a final experiment, 
Peng and co-workers successfully transferred nanopillar 
patterns first formed in PMMA onto an elastomeric film 
of poly(dimethylsiloxane) (PDMS) i.e. negative replica- 
tion of the original pattern. This demonstration outlined 
the ease with which potential patterns can be transferred 
into subsequent films for applications involving large area 
patterning. 



B. Motivation for this study 

In recent work |31j , we re-examined the prevailing hy- 
pothesis for pillar formation in nanofilms based on co- 
herent reflections of acoustic phonons in molten polymer 
nanofilms [261 I30j . Such a mechanism requires coherent 
phonon propagation of the order of the film thickness in 
an amorphous fiuid layer. A review of the literature has 
shown that acoustic phonon mean free paths of the order 
of 10-100 nm have only been measured in solid poly- 
mer nanofilms at frequencies of order 100 GHz and at 
temperatures -193°C < T < 27°C ^, far below the 
temperatures used in the experiments described above. 
Such long attenuation lengths, however, are highly un- 
likely in molten amorphous films far above Tg because of 
the degree of disorder present and the enhanced mobility 
of polymer chains at temperatures above Tg. 



Given that the free surface of thin liquid films is easily 
deformed by surface stresses }15j , we instead demonstrate 
in this work that nanopillar formations are caused by 
the nanoscale analogue of the long-wavelength Benard- 
Marangoni instability |33H36| , previously investigated for 
film thicknesses ranging from several hundred microns 
(70 "^ ho '^ 270 /im [311 inZj) to millimeters. In macro- 
scopically thicker films, film protrusions caused by ther- 
mocapillary flow are stabilized by capillary and gravi- 
tational forces, such that only gentle surface deflections 
are possible [37]. Onset of instability in such films re- 
quires that the inverse dynamic Bond number D^^^gj = 
7TATfiim/pff/io^ > 2/3(1 + F)-\ where p is the liq- 
uid density, 7^ = {dj/dTl, 7 is the liquid surface ten- 
sion, ATfiiin is the temperature drop across the liquid 
layer, F ~ {1 — k)/{D + k — 1) is an order one con- 
stant, D — do/ho, and k = fcair/fcmcit- Estimates corre- 
sponding to the experiments of Schaffer et al. and Peng 
et a/, indicate that D'^y'' > O(IO^) and G - 0(10""). 
These critical values lie far beyond the regime previously 
investigated by Vanhook and co-workers [36l [37] in which 
^orot - 0(10-1 - 1) and G - O(10-i - 10^). These 
estimates indicate that nanofilms dominated by thermo- 
capillary flow should always undergo instability. In what 
follows, we therefore propose an alternative mechanism 
to the acoustic phonon model to help explain the for- 
mation of elongated structures in liquid nanofilms sub- 
ject to a transverse thermal gradient. The analysis pre- 
sented here indicates that the experiments conducted 
by Schaffer et al. and Peng et al. provide a rare window 
into the dynamics of the less common long-wavelength 
Benard-Marangoni (BM) instability without interference 
from the better known short-wavelength (BM) instabil- 
ity, which gives rise to the beautiful cellular convection 
patterns often photographed. 

There is an additional feature worth emphasizing in 
Fig. [IJa). In the absence of a top wafer, a transverse 
thermal gradient can still be established in a film heated 
from below by natural or forced convection within the 
gas layer above the polymer melt. Since the Biot num- 
ber /3/io/fcinGit is linearly proportional to the polymer film 
thickness ho (where (3 is the heat transfer coefficient for 
natural convection), however, this number will be small. 
As a result, the thermal gradient within the viscous film 
will also be small and thermocapillary stresses at the in- 
terface may be easily stabilized by capillary forces. This 
is probably the reason why no fluid elongations were ob- 
served in the experiments of Chou et al. in which the 
polymer melt was heated in open air. Use of a top sub- 
strate maintained at a cooler temperature held in close 
proximity to the melt surface enforces a sizeable trans- 
verse thermal gradient which can be used to maximize 
and control thermocapillary flow. 

In this work we demonstrate that the predominance 
of thermocapillary forces along the free surface of molten 
nanofilms leads to a linearly unstable system which forms 
periodic protrusions no matter how small the applied 
thermal gradient in any liquid nanofilm, not just molten 



polymeric films. The analysis corresponds to a limit- 
ing case of Benard-Marangoni flow peculiar to viscous 
films of nanoscale dimensions such that hydrostatic forces 
are completely negligible and deformation amplitudes are 
small in comparison to the array pitch. Predictions of 
the pillar spacing from the linear analysis as a function of 
the substrate separation distance reveals good agreement 
with experiment. Deviations are likely due to overesti- 
mates in the reported values of ho for unannealed films, 
uncertainties in the measured values of do caused by the 
use of a tilted upper plate, and possible changes in wave- 
length caused by prolonged contact with the cooler sub- 
strate and film solidification prior to measurements of the 
array pitch. Finite element simulations of the full non- 
linear equation are also used to examine the array pitch 
and growth rates beyond the linear regime. Inspection of 
the Lyapunov free energy as a function of time confirms 
that in contrast to typical cellular instabilities in macro- 
scopically thick films, pillar-like elongations are energet- 
ically preferred in nanofilms. Provided there occurs no 
dewetting during film deformation, it is shown that fluid 
elongations continue to grow until contact with the cooler 
substrate is achieved. Identiflcation of the mechanism re- 
sponsible for this phenomenon may facilitate fabrication 
of extended arrays for nanoscale optical, photonic and 
biological applications. 



II. EVOLUTION OF MOLTEN NANOFILMS 

SUBJECT TO THE SLENDER GAP 

APPROXIMATION 

A. Films confined by parallel substrates 

The molten layer is modeled as an incompressible New- 
tonian fluid since the flow speeds and shear rates inher- 
ent in the experiments described are very small. Con- 
sistent with the slender gap approximation, all lateral 
dimensions are scaled by the pillar spacing distance L, 
while all vertical scales are normalized by the initial film 
thickness ho such that {X,Y) — {x/L,y/L), Z = z/ho, 
H{X,Y,t) — h{x,y,t)/ho and Do = do/ho- The pillar 
spacing L will later be identified with the wavelength of 
the maximally unstable mode, Amax, obtained from linear 
stability analysis. The conservation equations for mass 
and momentum within the thin liquid film are given by 
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Equation (pi) yields the scaling for the velocity compo- 
nents, namely if = {U,V,W) = {u/uc,v/uc,w/euc), 
where Uc represents the characteristic lateral speed set by 
thermocapillary flow. The corresponding Reynolds num- 
ber based on the initial film thickness is Re = pUcho/rj, 
where p and rj denote the polymer melt density and 
viscosity. In what follows, the polymer viscosity is as- 
sumed constant (i.e. a Newtonian fluid) and equal to 
7] — ?y(T'2) |38j . The non-dimensional Lagrangian or sub- 
stantial derivative is denoted by D /Dt = O/Ot +IJ-V 
where t = Uct/L. The overall (dimensionless) pressure 
in the fluid is given by 



P = eho{p + cj)) / {rjUc) 



(6) 



where p is the (dimensional) capillary pressure and 
represents contributions from hydrostatic pressure (i.e. 
(j) = g z where g is the gravitational constant) and dis- 
joining pressure (e.g. van der Waals forces). 

Within the slender gap approximation, e^ = 
{ho/ L)'^ ^ 1 and eRe — >■ 0; all terms on the left hand 
side of Eqs. (Isl) - ((sl) therefore vanish. In this limit, the 
pressure P within the thin film is independent of the ver- 
tical coordinate Z. Equations (pi) and (pi) can therefore 
be integrated with respect to Z , subject to the boundary 
conditions (BCs) at the liquid/solid and gas/liquid inter- 
face. Along the bottom substrate, it is assumed that the 
melt obeys the no-slip condition i.e. U\\ = {U, V) = 0. 
The dimensional stress jump across the air/melt inter- 
face |39j . which accounts for both normal and tangential 
stresses, is given by 



(Ta 



t)-n + V,7-7n(V, •n)==0. (7) 



Here, T = — (p + 0)1 + 27]E denotes the total bulk stress 
tensor, where I is the unit tensor and E the rate of strain 
tensor, n denotes the unit vector outwardly pointing from 
the melt interface, Vj represents the surface gradient op- 
erator [40! and 7 is the surface tension of the polymer 
melt in air. Since the viscosity and density of air are 
negligible in comparison to those of the melt, Tair = 0. 

Thermocapillary fiow within the melt leads to a non- 
vanishing shear stress Vs7 along the gas/liquid interface 
[39] . After a straightforward derivation, it can be shown 
within the slender gap approximation |35j that the tan- 
gential components of Eq. ([7| reduce to 

dU/dZ\z^„i^x,Y,r) - dr/dX (8) 



dV/dZ\z^HiX.Y.r)^dT/dY 



(9) 



where the surface gradient simplifies to Vg = Vy = 
{d/dX, d/dY). The variable F — e'-f/{r]Uc) represents the 
dimensionless surface tension. The gradients in surface 
tension arise directly from thermal gradients along the 
melt interface i.e. Vj|7 = [d'y / dT)V \\T . In dimensionless 
form, this relation is given by 
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where & = (T-Ti)/(T2-Ti), jt = \d"f/dT\, AT = T2- 
Ti > 0, and the Marangoni number Ma = e^T ^T / (rjUc) . 
In what follows, it is assumed that T2 — Ti > 0; further- 
more, for the liquid films of interest, the surface tension 
decreases linearly with increasing temperature T, which 
is reflected in the choice of the negative sign above. 

The in-plane velocity components are therefore given 
by: 
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Equation (111 represents a linear superposition of pres- 



sure driven flow caused by variations in interfacial curva- 
ture and hydrostatic forces, as described by Eq. ( 15 ), and 



shear driven flow induced by thermocapillary stresses. 
Substitution of Eq. (11) into Eq. ^ followed by in- 
tegration subject to the condition W(X, Y, Z = 0) = 
gives the vertical component of the velocity field, 
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The evolution equation for the moving interface can then 
be determined by integration of Eq. (l2| from < Z < 
H{X,Y,t) subject to W{X,Y,Z = 0) = and the 
kinematic boundary condition, W\z=h = DH/Dt = 
dH/dr + lJ\z=H ■ ^sH. The Leibnitz rule for differ- 
entiation gives 
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Substitution of Eq.(ll| leads to the evolution equation 



for the melt interface i/(X, y, r), namely 
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It is expected that the slender gap approximation re- 
mains valid throughout the growth process so long as 
{do/LY ^ 1, which holds for all the experiments de- 
scribed. 

Since the pressure in the film is independent of Z to 
order e^Re, one can determine its value by considering 
the normal stress balance at Z = H. The normal com- 
ponent of Eq. ([7| within the slender gap approximation 
yields the total pressure in the film to order e^: 



P = -Ca 



^^VniJ 



Ca'"^ BoH, 



(15) 



where Ca = rjUc/ije'^) and Bo = pgL"^ /^. Parameter 
estimates from the experiments of Schaffer et al. indi cat e 
that Ca^ is of the order of 10^ - 10^ [using Eq. ^] 
while Bo is of the order of 10^^ — 10^^. The hydrostatic 
contribution to the fluid pressure in Eq. (15) can there- 



fore be neglected altogether. The influence of disjoining 
pressure arising from van der Waals interactions in films 
ranging from 10 - 100 nm in thickness |35j is also ignored 



in this work. The flow induced by these molecular forces 
is weak in comparison to flow induced by thermocapillary 
stresses, which are of considerable magnitude in the ex- 
perimental systems of interest. While disjoining pressure 
effects can be included in straightforward fashion within 
P, they are not the primary mechanism for instability. 
Furthermore, there is yet no consensus in the literature 
on the appropriate analytic form of the disjoining pres- 
sure in cases where films are subject to large thermal 
gradients; most of the simplified forms available in the 
literature are only appropriate for isothermal systems. It 
is also assumed that any thermocapillary effects caused 
by solvent evaporation and subsequent cooling of the in- 
terface [41] can be neglected. This assumption requires 
that solvent evaporation be completed (either naturally 
or by film annealing) before the film is inserted into the 
experimental assembly. 

With these assumptions, the gradient of the Laplace 
pressure is given by 



ViiP = -Ca ^WrlH - e^VifiJViir. 



(16) 



The last term, which represents a correction to the 
Laplace pressure due to local variation in surface ten- 
sion, scales as e^ and can be safely ignored. The surface 
tension coefficient in the Laplace pressure only is there- 
fore set to the value 7 — j{T2). 

Determination of the interfacial stress conditions in 
Eqs. (Is]) and ^ requires knowledge of the thermal distri- 
bution along Z = H, which can be obtained from the en- 
ergy equations |35| pertaining to the confined air/liquid 
bilayer shown in Fig. [I] 



eRePr 



DO 
1)7 



d^e d^e 



dX^ dY' 






(17) 



Here, the Prandtl number Pr = v jot refers to the kine- 
matic viscosity v and thermal diffusivity a of the corre- 
sponding air or liquid melt layer. The Reynolds num- 
ber i?e, defined previously, is based on the correspond- 
ing layer thicknesses. Despite that Pr is of the order 
of 10^ — 10^ for the polymer melts of interest, the small 
gap approximation coupled with the vanishingly small 
value of Re (see Tables ll] and lll| ensures that the left 
hand side of Eq. (17 1 is completely negligible. In fact. 



the slender gap approximation is well satisfied in all the 
experiments described earlier since e^ ^ 1, ei?e ^ 1 
and eRePr <^ 1. The thermal analysis conveniently 
reduces to a one dimensional thermal conduction prob- 
lem for heat flow across an air/liquid bilayer subject to 
isothermal boundary conditions at Z — and Z — Do- 
The temperature distribution along the melt interface is 
therefore given by, e\z=H = (Do - H)/[Do + {n- l)H]. 
Substitution of this solution into Eq. ([lO| yields 



Viir = 



K Ma DoV\\H 



(18) 



Substitution of Eqn. (16 1 and (18 1 into Eq. (14) 



then yields the expression governing the motion of the 



TABLE I: Order of magnitude estimates for characteristic 
numbers used in the thermocapillary model extracted from 
the experiments of Schaffer et al. |22l [26] . Values of Pr for 
PS and PMMA at 170°C were obtained from Refs. [25] and 
|42| . The constant value of capillary number, Ca, results from 
the choice of thermocapillary velocity used to scale the flow 
speed, as discussed in the section following Eq. ( 30 1 



e 


10-^-10^2 


Re 


10"^**-10~^^ 


Pr 


10«-10^ 


Ma 


10'' -10^ 


Ca 


52.6 


Bo 


lo-'^-io-'' 



TABLE II: Literature values for air and polystyrene melt 
[M„ ^ 107 kg/mol and M^/M„ = 1.07 where n and w 
denotes number avg and weight avg] used in the analysis and 
numerical simulations. All parameter values quoted are for 
T = 170 °C, except for 7 and "fr, which were only available 
for T = 180 °C. For comparison, Schaffer et al. [22l |26] 
used polystyrene melts for which Mn ~ 108 kg/mol and 
M^/Mn = 1.03. 





Air 


PS 


p 


(kg/m^) 




0.829 [43] 


987(25] 


V 


(Pa- s) 




2.48 • lO-'^ [13] 


2.5 • 10* [12] 


k 


[W/(m°C)] 




0.036 [43] 


0.130 [25] 


a 


(mVs) 




4.25 • lO"'^ [43] 


6.45 • 10"* [25] 


1 


(10-=* N/m) 






31.53 [S] 


7t 


[10-^ N/(m 


=C)] 




0.0885 |44] 



air/liquid interface, namely 



dH 



K Do Ma H^ 

2[Do + {k-1)HI 



rVuH 



g3 

3Ca 



--VfiH 



(19) 

The characteristic scale for the lateral velocity, Uc, is 
set to the value established by thermocapillary flow. 



which can be obtained from Eq. ( 18 ) by letting the 



film thickness, slope and and interfacial stress be order 
one and equal to unity - i.e. H = 1, V||iJ = 1, and 
{dU/dZ)z=H = dV/dX = 1, such that 



Ma = 



{Do + n-lf 
nD„ 



(20) 



Since Ma = ejT ^T / (rjUc) , the scale for Uc becomes 
e K Do 7t at 



v{Do 



1) 



(21) 



The evolution of film disturbances governed by ther- 
mocapillary effects, as given by Eq. (19), is compared 
to evolution by acoustic phonon radiation pressure, as 
proposed by Schaffer et al. . While their derivation is 
also based on the slender gap approximation, the acous- 
tic phonon model neglects altogether any flow induced 
by tangential stresses due to interfacial thermal gradi- 
ents. Instead, the Laplace pressure, is counteracted by a 
radiation pressure due to phonon reflections which causes 
protrusions to grow. The overall fluid pressure in the AP 
model is therefore given by 

P= -Ca-^VlH -C^-^Q/[Do + {k-1)H], (22) 

where Q = 2Qkiii^AT / (upje^) and Q is the acoustic qual- 
ity factor described in the Introduction. Substitution of 
Eq. ([22J) into Eq. ([T4| (with V||r = since thermocap- 
illary effects play no role in the acoustic phonon model) 
yields the evolution equation proposed by Schaffer et al. :. 



dH 



g (1 - k)h^ 

3C^[Do + {k- l)Hf 



g3 

3Ua 



:V||iJ+-^=Vi^iJ =0. 

(23) 

Values for the thermophysical properties of air and PS 
are listed in Table |ll] Corresponding numbers for exper- 
iments with PMMA [53] are of similar magnitudes. 



Linear stability analysis of evolution equation 



Equations ( 19 ) and ( 23 ) can be further analyzed by lin- 



ear stability theory to provide an estimate of the fastest 
growing mode, the one most likely to be observed in ex- 
periment. Predictions of the corresponding wavelength 
are therefore expected to compare favorably with the 
pillar spacing measured in experiment if the proposed 
mechanism is correct. 



The behavior of Eq. ( 19 ) is examined in the limit where 
an initially flat and uniform film of thickness H = 1 
(i.e. base state) is subject to an infinitesimal periodic 
perturbation of amplitude SHo ^ 1 and wave number 
A II where \K\\\ = K = 2TrL/\. Solutions of the form 
H{X,Y,t) ^1 + SHo exp[/3(if)r]exp[i^|| -^n] are sub- 
stituted into Eq. (19), where X\\ = iX,Y), and all 
quadratic or higher order terms are neglected. The re- 
sulting expression for the growth rate is 



/3W = 



KD„AIa 



"' ] K-'. 



2{Do + n-lY iCaJ 



(24) 



Disturbances for which P{K) — neither grow nor decay. 
This condition establishes the criterion for marginal (M) 
stability where the corresponding wave number, Ky^^ for 
the thermocapillary model, is given by 



K, 



TC 

M 



uDoMa Ca 



2 {Do 



ly 



(25) 



Note that in the absence of any stabihzing hydrostatic 
terms as is the case with nanofilnis, there always exists a 
band of wavenumbers < K < Kjjp for which the film 
is linearly unstable, no matter how small the value of the 
imposed temperature gradient. This stands in sharp con- 
trast to the thermocapillary instability in much thicker 
films IMl ^ for which /Cm = (3 k DoW^C^/[2{Do + 
K — 1)^] — BoY^^. For thicker films, there exists a critical 
Marangoni number for onset of instability: 



Mar, 



2 Bo (L>o + K - 1)2 

3 Ca kDo 



(26) 



This criterion 
of the 
7TATfii„/(pff/i2) 



is commonly expressed in terms 

dyn _ 

onset 

where ATk 



inverse dynamic Bond number D, 



> 2/3(1 + F)- 



Cfilm 



represents the temperature drop across the liquid layer 
[and not the temperature drop across the bilayer as 
defined in Eq. ([26])] and F = {I - k)/{Do + k - 1) is a 
constant of order one [37] • The regime investigated by 
vanHook et al. for films of the order of several hundred 
microns corresponds to values of -Donsot ^^ ^^^ range 
10~^ — 1. By contrast, representative values for IJ'^y" 
in the experiments of Schaffer et al. and Peng et al. are 
of the order of 10^. Nanofilms subject to a transverse 
thermal gradient are therefore always linear unstable 
irrespective of the magnitude of AT. 

The fastest growing wave number is determined from 
the extremum of f3{K) in Eq. (24), with the result that 
-^max = -^M^/>/2 = 27rL/A„iax. In dimensional units, 
the wavelength of the most unstable mode is given by 



\TC 



2TTho\ 



A-fho 



3 ndojT AT 



do_ 
hr, 



K-1 



(27) 



This expression provides an estimate of the average spac- 
ing between protrusions undergoing growth by thermo- 
capillary flow. For the nanofilm experiments described 
earher, ho ~ O (100 nm), ho < do < 8 ho and AT w 
10 — 50°C. This leads to predictions of the pillar spacings 
ranging from about 2-20 /im. (More detailed comparison 
to experiments will be discussed in Section III.) Accord- 
ing to Eq. (27), the characteristic lateral spacing be- 



tween nanopillars is determined by the initial film thick- 
ness, ho, as well as the gap ratio Do = do/ho, the ratio 
of the surface tension to the maximum change in surface 
tension, 7/(72^ AT), and the ratio of thermal conductivi- 
ties K — fcair/fcmcit- For cascs in which the geometry and 
material properties are held fixed, a larger thermal gra- 
dient produces more closely spaced pillars. Reversal of 
the thermal gradient such that T2 < Ti should lead to 
linearly stable films. 

Figure W[a) represents solutions to Eq. (27) for a 



polystyrene film at T2 = 170°C with AT = 43 °C. 
Smaller gap ratios Do lead to smaller values of pillar spac- 
ing since the film is subject to a larger effective thermal 
gradient. Figure [2jb) highlights the dependence of A™^ 
on the initial film thickness ho for various gap widths do 



and AT = 43° C. As evident, the prediction for Aj^^x de- 
pends sensitively on ho, especially for the smallest values 
of hn- 
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FIG. 2: Solutions of Eq. ([27| for k = 0.277 and AT = 43°C. 
Curves show a sharp decrease in A^Jx for the smaller values 
of ho. 



The linear stability analysis of Eq. ( 23 ) yields a predic- 



tion for the fastest growing wavelength for the acoustic 
phonon model, namely Eq. (11]). The ratio of dominant 
wavelengths corresponding to the two proposed mecha- 
nisms is given by 



A^ ^ / 4Qfcniclt(l-K) 
A^ax V 'i Up -IT Do 



(28) 



Future experiments conducted with parallel substrates 
for a wider range of Do should help identify the operating 
mechanism leading to pillar formation. 



The characteristic velocity defined earlier in Eq. (21 1, 



which sets the scale for the lateral flow speed based on 
thermocapillary stress, can be re-expressed in terms of 
the length scale Aj^^x obtained from linear stability anal- 
ysis: 



XA^f 



ATC 



7 
77 



(4^) 



(29) 



Here, the lateral scale L used to define the slender gap pa- 
rameter, e = ho/L, is identified with Aj^^x- Similarly, the 
characteristic timescale based on thermocapillary flow is 
given by 



aTC 
max 



3/i-o 

(47r)2 



(30) 



Estimates from the experiments of SchafFer et al. indicate 
that Uc is of the order of 10~^ — 10^ nm/s and ic ~ 
0(10^^ — 10'^ hrs). If the thermocapillary flow speed u^ 
given by Eq. ( 29 ) is used to define the capillary num- 
ber, then Ca = (47r)'^/3, a fixed constant. Replacing the 
capillary number by this numerical value and substitut- 
ing the expression for the Marangoni number given by 
Eq. ( 20 ) into the interface equation Eq. ( 19 ) yields the 



following form of the evolution equation: 



dH 



Dq+k-I 
Do + {k-1)H 



n 2 



H' 



-V\\H 



H-' 



rViH). =0. 



2 ■ "" ■ (47r)2 

(31) 
For thicker films, hydrostatic forces can be re- 
incorporated into this expression by including the term 
—Bo H^V\\H/{A-kY in the curly brackets. During the 
early stages of film deformation when H and VyiJ are 
order one, the relative magnitude of terms in Eq. (31) 



reveals the basis for pillar formation. The ratio of ther- 
mocapillary to capillary flux scales as Stt^, while the 
ratio of thermocapillary to gravitational flux scales as 
Stt^/Bo w 10^ — 10^. These estimates reveal that ther- 
mocapillary forces overcome the stabilizing effect of cap- 
illary and gravitational forces even at early times. In 
section IV. B, it is shown that thermocapillary forces pre- 
vail even more strongly at late times for parameter values 
pertinent to the nanofilm experiments. A similar com- 
parison can be made using the parameter values in the 
experiments of VanHook et al. j45) with thicker films 
(70 ^ ho '^ 270 fim) and much smaller transverse ther- 
mal gradients (180 < Ta - Ti/do < 500 °C/cm). While 
the thermocapillary to capillary flux ratio remains at Stt^, 
the thermocapillary to gravitational flux ratio decreases 
to 10~^, eight to nine orders of magnitude smaller than 
the ratio in the nanofllm experiments of Schaffer et al. 
and Peng et al. . While gravitational forces effectively re- 
press the growth of pillars in macroscopically thick fllms, 
this order of magnitude analysis confirms that hydro- 
static forces are ineffective in repressing the growth of 
elongations in nanoscale films. 



Integration of the full nonlinear Eq. ( 19 ) can be used 



to compute a lower bound on the time interval, itop, 
required for nanopillars to contact the cooler substrate 
within the approximation of a constant film viscosity 
[55] . It will be shown in Section IV. A that estimates ob- 
tained from the growth rate of the most unstable mode, 
/3(^max)j are in fairly good agreement with the estimates 
obtained from numerical solutions of Eq. ( 19 1 for the pa- 
rameter range of interest. Substitution of Eq. (201 and 
Ca = (47r)2/3 into Eq. ^ and Eq. ^ yields the sim- 



plified expression for the growth rate: 



K' 



(32) 



The wave number corresponding to marginal stability is 
therefore KJ[^ = A^:/^/2. Since K'I^^ = Klp/^/2 = 27r, 
the growth rate for the fastest growing mode simply re- 



duces to /3, 



TC 



Setting 5Ho exp[/3( JCn 



leads to the expression Ttop 



ln[(D„ 



l)/SHo/7r\ 



which in dimensional units corresponds to ttop = 

(3 77/io/7)[A^c^/(27r/io)]'*ln[(i:'o - l)/SHo]- Substitution 
of Eq. (27 1 into this expression then gives 

16TJ-fho{Do + K-l)* 



^top ~^ 



3 {k Do it Ar)2 



In 



Do -I 
SHo 



(33) 



Estimates of ttop for the nanopillar experiments range 
from about tens of minutes tjjjens of hours for the 
largest gap spacings used and SHo = 10~^. Low molec- 
ular weight polymers with much smaller viscosities re- 
quire proportionally less time to contact the cooler top 
substrate. Studies of this sort are useful in determining 
when to remove the thermal gradient in order to form 
nanopillars of specified height. 



2. Lyapunov free energy for evolving interface 

Hydrodynamic systems subject to interfacial insta- 
bility sometimes exhibit steady states as observed in 
Rayleigh-Benard or Benard-Marangoni cellular convec- 
tion. Within the context of the experiments described, 
this would require pillar formations which once formed, 
neither grow nor decay, representing a fixed spatial con- 
figuration while the melt continues to undergo surface 
and interior flow. To examine this possibility, one can in- 
vestigate the temporal behavior of the Lyapunov free en- 
ergy associated with the evolving interface, as previously 
implemented in Re fs. [iGlll^. This approach is based on 
the analysis of interface problems using the well known 
form of the Cahn-Hilliard free energy for systems with 
spatial variation in an intensive scalar variable like com- 
position or density [IS]. The Cahn-Hilliard equation has 
been successfully used to explore the evolution of moving 
interfaces in binary systems undergoing phase separation. 
This approach, which involves monitoring the free energy 
associated with the entire film undergoing deformation, 
provides a more accurate assessment of possible steady 
state configurations than simple considerations based on 
Eq.(|3l]) in the limit dH/dr -^ 0. 



In the Appendix, it is shown that the free energy 
corresponding to the nanofilms of interest is given by 
5 = / £ dXdY, where 



£=(V||iJ) 



SnMa Ca 

D'o 



iJln 



H 



I + XH 



+ln(l + x) 



(34) 



and X = ('^ ^ l)/-Do. Numerical solutions of Eq. (A 
13 1 for large and small values of the gap ratio, Dg, are 
discussed in Section IV. B. 
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B. Films confined by non- parallel substrates 

The analysis presented in Section II. A describes the 
evolution of a fluid bilayer interface confined by two flat 
and parallel substrates separated by a distance Do — 
do/ho- As described in Section I. A. 2, however, Schaffer 
and co-workers purposely used in all their experiments a 
tilted plate geometry in which the top wafer was inclined 
with respect to the bottom one by about 1/Lim/lcm, cor- 
responding to an inclination angle ip of about 0.0057°. 
The evolution equation can be modified to account for 
two flat substrates with relative tilt. When the cooler 
substrate is tilted away from the horizontal by a con- 
stant angle ip, the local value of the plate separation will 
depend on iX,Y) such that D{X,Y) = d{x,y)/ho. This 
modification alters the film surface temperature, Q\z=h, 
as well as the surface thermal gradient, WQ\z=h, which 
in turn alters the interfacial thermocapillary stress, ViiF. 
Accordingly, 



e\z=H^{D-H)/[D + {K~l)H] 



(35) 



Vi,ei 



Z=H 



[D+{k-1)H]'- 



{HWiiD-DWi\H) 



(36) 



and 



Vnri 



Z=H 



-MaV\\e\z=H 



KMa 



[D+{k-1)H]^ 



(37) 



Here, D(A ||) = Do+t&n{ip)T \\ '-^ h '^^'^re Dg represents 

the gap ratio at Ay = (later identified with the mid- 
point of the computational domain). The quantities ^ 
and tan^ = tamp/e represent variables rescaled accord- 
ing to the slender gap approximation. In the numerical 
solutions discussed in Section IV. C. 2, the tilt of the upper 
substrate is defined by the unit vector Ty = (1, 1)/a/2. 
Substitution of Eq. 0f\ into Eq. ([Tel) and Eq. (fwl) leads 



to the modified evolution equation 



dH 



Vii -Q 



tilt 



, 



(38) 



where 



^tilt 



KMaH'^{DWuH-Htmi{Tpyf») H^ 



2[D + [n - X)HY 



3Ca 



S/tH. 



(39) 
A linear stability analysis of Eq. ( 38 ) (not shown here) 



confirms that the pattern wavelength in Eq. (27) remains 



unaffected by the small tilt angle used in the experiments 
of Schaffer et al. . More generally, Eq. ( 27 1 remains valid 

so long as | tan(^)| < 0{5Ho). 



III. NANOPILLAR SPACINGS: COMPARISON 
BETWEEN EXPERIMENT AND THEORY 



Shown in Fig. l3Fa) is a direct comparison of Eq. (27) 
with the experimental data of Schaffer et al. [22l [26] The 
solid lines denote the predictions of the thermocapillary 
model with no adjustable parameter values using the ma- 
terial properties listed in Table [Til the symbols denote the 
experimental data. 
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FIG. 3: Direct comparison of theoretical estimates for dom- 
inant instability wavelength, Amax for both thermocapillary 
(TC) and acoustic phonon (AP) models with experimental 
measurements from Schaffer et al. [221 1261 1271 I30j as func- 
tion of increasing wafer separation distance do. (a) Plots of 
Eq. (271 for thermocapillary model with no adjustable pa- 
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rameters for different experiments labeled A-D using material 
constants listed in Table p] (b) Plots showing least squares 
fits to the TC and AP models. The TC model was fitted to 
the form given by Eq. (271, namely AJ^^x = Cwfd'o-^C^I \fd'o. 
AP model was fitted to the function given by Eq. ^\ with 
itp — 1850 m/s and Q — 6.2 as fitting parameters. Table [lll[ 
lists the fitting coefficients, C\ and C2, obtained for the TC 
model. 

While the overall functional behavior of A^^x with do 
is in good agreement with experiment, the model system- 
atically overestimates the pillar spacings, in some cases 
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TABLE III: Coefficients obtained from a least squares fit to 
the data of Schaffer et al. [221 126] . Experimental data were 
fit to the function A^ax = Ci do^''^ + C2 do~^'^ given by 



Eq. (27 1, where Ani^x is reported in microns and ho and do 
in nm. The constants C^'* denote the values obtained for 
the least squares fits shown in Fig. Isl The constants C"^ 
represent the predictions of the TC model given by Eq. ( 27 1 
with no adjustable parameters using the material properties 
listed in Table |ll| The percentage errors are defined by 

Pit ^TC 



{c- 



C''')/C 
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D 




ho (nm) 


80 
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100 
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AT (°C) 


43 
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cr 


[(10^ Mm)°-5] 
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0.38 
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cr 
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% Error Ci 




-0.53 


16 


1.2 


39 


% Error C2 




-69 


-21 


-66 


42 



by as much as 40%. This is especially evident in exper- 
imental run B for which ho = 96 nm and AT = 11 °C. 
Before discussing these discrepancies in detail, it is use- 
ful to examine a least-squares fit of the data to the func- 
tion A;^C^ = Ci do^/^ -I- C2 do~^^'^ given by Eq. ^, as 
shown in Fig. ^h). Listed in Table [!!!| 



IS a comparison 



of the analytic expressions for the two constants, namely 
Ci^c = 27r[4 ho 7/(3 K -fT Ar]i/2 and Cj'=^ = Ci{K-l)ho, 
along with the results for the fitting constants denoted 
by C^'*. 



In general, the agreement between the TC model and 
experiment improves for larger values of AT. However, 
given that the least squares fit captures the experimental 
trend with increasing values of ho and do so well, it is 
worth considering what experimental challenges might 
affect the reported measurements. For completeness, we 
include in Fig. [3[b) two additional dashed lines for runs 
B and C, which represent a least squares fit of the data 
to Eq. (IT]) with Q — 6.2 and Up = 1850 m/s, the same 
fitting constants reported by Schaffer et al. [351 [IS] 



A. Possible causes of discrepancy between theory 
and experiment 

There are several experimental challenges in perform- 
ing the experiments on nanopillar formation. Perhaps 



the most important is that all experiments to date have 
used silicon wafers to confine the polymer films. These 
opaque substrates prevent observation of the instability 
in-situ. In fact, measurements of the pillar spacings were 
normally obtained long after the pillars had contacted 
the cooler wafer. The pillar amplitudes were by then size- 
able, possibly violating the assumptions of linear stability 
analysis. Furthermore, the warmer nanopillars had sus- 
tained prolonged contact with a cooler substrate leading 
to possible reorganization of fluid due to thermocapillary 
or other packing effects along the underside of the top 
wafer. Measurements taken once the pillars had solidi- 
fied and the top wafer was removed may therefore differ 
from the predictions of linear stability theory. In many of 
the experiments described earlier, measurements of the 
spacing between fluid elongations included not only pillar 
arrays, but lamellar, spiral and other periodic structures 
since these were more commonly obtained. An additional 
complication is that a typical molten nanofllm is not com- 
pletely smooth and flat due to the presence of contam- 
inant particles and pinholes caused by dewetting. Any 
small fluid elevations caused by these nucleation points 
are prone to rapid growth when subject to a thermal 
gradient. Structures arising from such initial conditions, 
however, correspond more to disturbances of finite am- 
plitude and not infinitesimal amplitudes as assumed by 
the linear analysis. 

As evident from the curves in Fig. [2] the param- 
eters ho and do strongly affect the predicted values of 
Aj^^x- The sharp drop in Aj^^x becomes even more pro- 
nounced for smaller values of AT [3T]. Validation of 
either mechanism proposed therefore requires accurate 
measurements of the film thickness. It appears that the 
films used by Schaffer et al. ^^^^ and Peng et 
al. |23| were not annealed prior to insertion in the exper- 
imental setup. Spun cast polymer films tend to retain a 
significant amount of solvent [IHl 120] , which is normally 
expelled by film annealing in vacuum at elevated temper- 
atures for several hours. (Annealing has the additional 
advantage of healing pin holes that sometimes form dur- 
ing spin coating.) Significant film shrinkage typically ac- 
companies this process due to solvent evaporation. The 
degree of film shrinkage depends on the ambient vapor 
pressure as well as the time and temperature of the bake. 
It is therefore likely that the values of ho reported in the 
literature represent overestimates of the initial film thick- 
ness ho- Smaller values of ho lead to smaller predictions 
for the pillar spacing, in closer agreement with experi- 
ment. 

The distance between pillars in experiment was typ- 
ically obtained by direct measurement from optical mi- 
crographs. In future experiments, it would be prefer- 
able to Fourier analyze the patterns obtained by an EFT 
(Fast Fourier Transform) analysis. This analysis may 
reveal not only the dominant wave number but harmon- 
ics that develop due to the growth of smaller pillars in 
between two larger neighboring ones. Such an analysis, 
however, requires a fair number of protrusions for sta- 
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TABLE IV: Typical values of sensitivity coefficients S^. /A 
resulting from Eq. 1 27 1 



TC 
max 



where £,i =7, -yr, AT, k, Do or ho. 



i'7/Amax 


0.5 


S-yy/Aniax) SAT/Aniax 


-0.5 




-0.5 to -0.25 


SDo /Amax 


0.6 to 1.1 


'-'^o/'^max 


to 0.4 



tistically meaningful results. It may have been the case 
with the tilted plate geometry, that the smaller domains 
corresponding to each distinct value of Do forbade use of 
this technique. 

We conducted FFT analyses of nanopillar arrays 
published in the literature |551 [23 HZj and were sur- 
prised to find a very wide distribution in pillar spac- 
ings within even a single experiment. Often there ap- 
peared not a single dominant wavelength but several 
competing wavelengths. This finding prompted a sen- 



sitivity analysis of Eq. (271 to better understand which 



variables most strongly affect the uncertainty in mea- 
surements of AJ^^x — ^nia.xi^i)j ^ defined by Uxtc = 
^y'^^iS^. A^j/^j)2. Here, the relative sensitivity coef- 
ficients are given by 5*^. = ^i dX]^^/d$i where £_i = 
(7, 7t, at, k. Do, ho). This analysis demonstrates that 
Sj == -S'-yy = -Sat = 1/2 A™^, S,, = [n/iDo + k- 
1) - 1/2]ATC^, Sn^ - [Do/{Do + k - 1) - 1/2]ATC^ and 
Sh^ = [1/2 + {k- 1)1 {Do + k- 1)]A^2^. Typical values 
for these sensitivity coefficients for the parameter values 
corresponding to the experiments of Schaffer et al. are 
summarized in Table IIVI These values indicate that the 
gap ratio. Do = do/ho, the initial film thickness, ho, and 
the polymer surface tension, 7, most significantly influ- 
ence the degree of uncertainty in measurements of A 



TC 
max* 



IV. NUMERICAL SIMULATION OF THIN 

FILM EQUATION: LINEAR AND NON-LINEAR 

REGIMES 



forced along the domain edges (except for the simulations 
using tilted substrates). A quadrilateral mesh consisting 
of 200 X 200 elements was applied for the coarse (non- 
extended) discretization, leading to an extended system 
of equations with about 5 • 10^ degrees of freedom. An 
imphcit Newton iteration scheme was used to advance 
the position of the film interface in time; the linear sys- 
tem of equations for each iteration was solved using the 
iterative solver GMRES (Generalized Minimal Residual 
Method). All simulations were conducted on HP Pro- 
Liant DL360 G4p workstations equipped with dual Intel 
Xeon 3.0 GHz processors running CentOS 4.6. The typ- 
ical growth of a nanopillar spanning two substrates (i.e. 
T ~ Ttop) required approximately 5 — 6 hrs of CPU time, 
corresponding to about 900-1000 integration steps. Nu- 
merical convergence tests were conducted by evaluating 
the local dimensionless film height at A^ = 400 inter- 
polation points within the square domain. These tests 
confirmed that both the average difference, AiJavg = 
Eti \H2{X,Y,T,,^) - Hi{X,Y,T,,^)\/N, as well as the 
maximum difference, AiJmax = max |_ff2(^, ^, '''top) — 
Hi{X,Y,Ttap)\/N , in film height at the end of a run 
(i.e. T ~ Ttop) were less than 10~^ when decreasing 
the grid size or integration time step. Here, Hi denotes 
the coarser measurement and H2 the refined one. Fur- 
ther tests revealed that the film volume was conserved 
during each run to a value l^V/V = \ ^[H{X,Y,Ttop) — 
H{X,Y,T^Q)] dXdr|/(AXAy) < 10^1". 

In all simulations conducted, the thickness of the ini- 
tial flat film was modulated by a very small amount of 
white noise such that H{X, Y, r = 0) = 1 + ^ IH, where 
d\ denotes a random number between —1 and +1. The 
amplitude of the white noise was set to 1^ = 0(10^^). Ac- 
cording to Eq. ( 33 ), larger values of ^ will lead to shorter 
contact times in proportion to —\ti{5Ho). In order to 
facilitate detailed comparison between runs for different 
choices of experimental parameters, the random number 
algorithm was reset before each run so as to generate 
an identical white noise distribution. Initialization with 
white noise was preferable to initialization by a sinusoidal 
function, as is common, in order not to bias the system 
toward a preferred wavelength too early in the pillar for- 
mation process. 



To investigate the extent of non-linear effects on the 
growth of nanopillars, we also conducted 3D finite ele- 
ment simulations of Eq. ( 19 1 using a commercial soft- 



ware package [511. Material properties corresponding to 
molten polystyrene (PS) were used in these numerical 
studies (see Table Ih]). The computational domain corre- 
sponded to a square of size AX x AF = OAj^^x ^ 6A™x 
(according to Eq. (27)) where spatial discretization was 
obtained via second order Lagrangian shape functions. 
This choice in domain size and discretization order re- 
fiects a compromise between available computational re- 
sources and generation of a sufficient number of peaks 
for FFT analysis. Periodic boundary conditions were en- 



A. Films confined by parallel wafers 

FFTs of the in-plane images obtained from the numer- 
ical solution of Eq. (19 1 were used to extract values of 
the dominant wavelength, A^™"'(t), at each instant in 
time. This numerical value was co mpa red to the theoret- 
ical prediction A^J^^^ given by Eq. (27). Shown in Fig. H 
are results of these simulations. 

The FFTs were computed by sampling 200 x 200 points 
within the computational domain for each value of r; ap- 
proximately 140 instances in time were so evaluated. The 
legend in each plot represents the variables held fixed dur- 
ing the simulation; the table entries specify the theoreti- 
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FIG. 4: Direct comparison of Ajnax from Eq. (271 with the 
instantaneous wavelength, A™x'i extracted from FFT analy- 
sis of numerical solutions of the evolving film thickness from 
Eq. ( |19[ ) with increasing dimensionless time r. Symbol rtop 
represents time of contact of the fastest growing pillar with 
cooler substrate. Times Tq = 0, ti = 0.06 and rtop ~ 1-23 
shown in the top panel refer to time stamps of the snapshot 
images shown next in Fig. [s] (a) Variation of the wavelength 
ratio with increasing temperature difference AT. (b) Varia- 
tion of the wavelength ratio with increasing gap separation 
distance do- (c) Variation of the wavelength ratio with in- 
creasing values of initial film thickness ho. 



cal values of Aj^^x corresponding to the chosen parameter 
set. For convenience, the factor used in converting r to 
real time t is also listed. The times Tq ~ Q, ti = 0.06 
and Ttop = 1-23 shown in Fig. W[b,) denote the three in- 
stances in time for for which the FFTs shown in Fig. [s] 
were computed. The variable rtop denotes the time at 
which the fastest growing nanopillar in a particular run 
made contact with the cooler substrate, at which point 
the simulation is terminated. The times rtop = I, II or III 
indicate this contact time for the parameters values des- 
ignated by (I), (II) or (III). 

As evident, the overall deviation of A^™"'(r) from Aj^^x 
is rather small regardless of the parameter range used. In 
all cases, this ratio rapidly approaches unity as r — )• 1. 
The only discernible difference is that the fastest grow- 
ing peaks require a longer time to contact the cooler 
substrate for larger values of the relative gap spacing 
Dq = do /ho, as expected. The very short lived but large 
initial transients are caused by initialization with white 
noise; shortly following r = 0, there exist disturbances of 
all wavelengths. Those contributions with wave number 
larger than the cut-off wave number K^ become rapidly 
damped. The ratio A™"x'(''')/'^inax then drops sharply to 
a value close to one as the maximally unstable distur- 
bance is established. The approach to unity from below 
rather than above is due to the asymmetry in the disper- 
sion curve (3{K) for which there exists a broader band of 
unstable wave numbers below i^max than above. 

Additional simulations (not shown for brevity) reveal 
that A^™"'/AJ^^x — >■ 1 by r = 1 irrespective of the spe- 
cific initialization function used i.e. white noise or a sim- 
ple sinusoidal function. Initialization by a double cosine 
wave in {X, Y) with wavelength Ac = X]^^/\/2, for ex- 
ample, produced the same long time behavior shown so 
long as the amplitude of the disturbance function satis- 
fied SHo < 1. 

Images of the evolving film thickness, H{X,Y,t) — 1, 
as seen from above, the corresponding Fourier trans- 
form (insets), and cross-sectional views along the mirror 
planes X = and Y — are shown in Fig. [s] at times 
r = 0, 0.06 and 1.23. The relevant parameters values are 
ho = 100 nm, do = 285 nm and AT = 46 °C, which rep- 
resent case II) in Fig. H. The arrow shown in the FFT 
with unit length denotes the magnitude i^max- ^^ evi- 
dent from the images in Figs.lsla) and (b), although the 
disturbance heights of order 10~^ do not increase sub- 
stantially from r = to 0.06, an increasingly regular 
hexagonal pattern is already visible, both in the Fourier 
transform as well the cross sectional views. Figure pfc) 
depicts the in-plane symmetry in a fully evolved film, just 
as the fastest growing peak contacts the cooler substrate. 
Here, the pillar amplitudes have increased substantially 
in comparison to their initial values. By this time, the 
Fourier transform of the emerging pattern has evolved 
from a wide band into a narrow ring with distinct six-fold 
symmetry and mean radius K^^^. Values of the (dimen- 
sionless) interfacial shear stress, Tx = dT/dX, along the 
axis Y = for r = 1.23 are shown in the bottom right 
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image. As expected from symmetry, the local extrema in 
film thickness along Y = (solid black curve) occur at 
the locations of vanishing shear stress i.e. Tx ~ 0. The 
largest values of |rx| tend to occur near the maxima and 
minima in film thickness. 
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FIG. 5: View from above of numerical solutions of the film 
thickness, H{X,Y,t) — f.O from Eq. (f9l at three instants 
in time: (a) tq — (origin of time), (b) ri — 0.06 and 
Ttop = 1.23. Simulation parameters values are ho= 100 nm, 
do=285 nm and AT — 46°C. Evolution of the correspond- 
ing dominant wavelength is depicted by case (II) in Fig.WTa). 
Left panel depicts amplitude H{X, Y, r) — 1.0 (white = eleva- 
tions, black = depressions); right panel depicts cross sectional 
views along axes X = and Y — 0. Inset images show the 
instantaneous 2D Fourier transform of the corresponding film 
thickness. Unit arrows denote magnitude of most unstable 
wave number, K^-^^^/\/^ = 2n, derived from linear stability 
theory [see discussion following Eq. (32l]. Values of the di- 
mensionless interfacial shear stress, Fx = dT/dX , along the 
axis F = are shown in the bottom right image. 

Shown in Fig. [g] is the growth rate ratio, /3n™x^V/5m£y 
for the parameter values labeled (b) in Fig. Hand Fig. [5] 
This ratio was computed for each of the six most rapidly 
growing peaks according to 
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FIG. 6: Direct comparison of instability growth rate, /3n 

TT^ (as discussed in Sec tion II.A.l), with instantaneous growth 

rate, /3^^', from Eq. (401, with increasing time r. Different 



curves shown correspond to growth rates of six fastest peaks 
for parameter values ho = 100 nm, do = 285 nm and AT = 
46° C. Inset image depicts film shape for H{X, Y, Ttop ~ 1.23). 



As in the solutions shown in Fig. |4J here too the nu- 
merical results are initially influenced by the white noise 
disturbance spectrum. Each of the six fastest growing 
peaks behaves somewhat differently at the earliest times 
depending on what is the local value of the disturbance 
height. However, the growth rates collapse rapidly by 
about T = 0.4, after which the average growth rate slowly 
increases toward the prediction of linear stability theory, 
which is established by about t — 1.0. Beyond this time, 
the solutions reveal rapid growth and an increasing de- 
parture from the predictions of linear stability theory as 
nonlinear effects contribute to the evolving pattern. Be- 
yond T « 1.0, the growing nanopillars are within reach 
of the cooler substrate. The instances marked tq, ti and 
Ttop represent exactly those times indicated in Fig. Ba) 
and Fig. [Sj 



B. Numerical simulations of Lyapunov free energy 



Numerical solutions of Eq. ( 19 ) confirm that nonlinear 



effects for the parameter sets examined become signficant 
only when fluid elongations come into close proximity 
with the cooler substrate. As evident in Fig. |6j the 
elongation rate then exceeds exponential growth. In this 
regime, the nanopillars have grown a distance large in 
comparison to the initial film disturbance heights and 
the nonlinear terms in Eq. ( 19 1 strongly influence the 



flow. To explore the energetics of formation beyond the 
linear regime, we investigated the temporal behavior of 
the Lyapunov free energy given by Eq. ( |34[ ). Shown in 
Fig. n\ are solutions of the free energy ^ — J £, dXdY 
for a polystyrene nanofilm with ho = 100 nm and AT = 
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46tl for two different wafer separation distances, do = 
285 nm and 800 nm. The termination points represent 
Ttop . The individual contributions to the total free energy 
(denoted by "Sum") from capillary and thermocapillary 
terms feature several important points. 
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FIG. 7: Numerical solutions of Lyapunov free energy, §(''"), 
for an initial flat film of thickness ho = 100 nm and temper- 
ature difference AT = 46° C subject to two wafer separation 
distances: (a) do = 285 nm and (b) 800 nm. 

For T < 1, the film experiences small deformations such 
that the opposing capillary and thermocapillary contri- 
butions are also small, neither significantly enhancing nor 
depicting energy from the evolving film. Magnified views 
of the curves (not shown) confirm a small but monoti- 
cally decreasing value of the free energy due to the still 
dominant influence of thermocapillary stresses. This pe- 
riod of growth corresponds to the linear regime described 
by linear stability analysis. Strong departure from this 
behavior occurs for t > 1 when nonlinear effects begin 
to dominate. In this regime, the time (or distance) re- 
maining for fiuid contact with the top wafer is small and 
the energetics of pillar formation strongly affected by the 
presence of the cooler target. For the smaller gap separa- 
tion distance (do = 285 nm) shown FigJ7](a), thermocap- 
illary effects dominate capillary effects as the nanopillars 



grow ever more rapidly toward the cooler target. There 
remains sufficient fluid in the residual film to continue 
feeding the growth of nanopillars such that the system 
continuously lowers its overall free energy by transport- 
ing fluid toward the cooler substrate. Unlike the equilib- 
rium cellular convective patterns observed with Rayleigh- 
Benard or Benard-Marangoni instabilities, this nanofilm 
instability is non-saturating and the free energy contin- 
ues to decrease until the fluid makes contact with the 
cooler target. 

The results shown in Fig. [71(b) for the larger gap sep- 
aration distance do = 800 nm reveal different behavior. 
Since the top substrate is positioned further away, the 
initial thermal gradient is smaller and the films require 
correspondingly longer times to develop substantial fluid 
elongations. The linear to nonlinear transition is ob- 
served to occur at slightly later times, r w 1.2. The indi- 
vidual contributions to the free energy are still clearly 
distinguishable but eventually asymptote. The larger 
wafer separation distance allows for longer growth peri- 
ods, which causes significant film depletion near the base 
of nanopillars. Fluid transfer needed to grow the elonga- 
tions is impeded, eventually halting their growth. Fluid 
already contained within the nanopillars continues to un- 
dergo a circulatory flow pattern, rising upwards near the 
surface due to thermocapillary stresses and falling down- 
wards near the interior due to capillary stresses. How- 
ever, fluid transfer from the initial deposited film slows 
considerably and can be halted completely if the deple- 
tion effect causes dryout. 

In summary, the Lyapunov analysis demonstrates why 
there is no steady state configuration in nanofilms ex- 
cept in cases where film depletion leads to pillar isola- 
tion. This limit can be achieved by placing the secondary 
plate sufficiently far from the initial deposited film. In 
this case, nanopillars that form will continue to undergo 
surface and interior flow but they cannot grow substan- 
tially in height due to a limitation in the available fluid 
mass needed to feed continued growth and elongation. 



C. Influence of relative gap spacing and substrate 
tilt on symmetry of evolving films 

1. Effect of larger gap spacing 

It is interesting to explore further the nonlinear behav- 
ior shown in Fig. p] for times t > 1 by examining images 
of the evolved films. The nonlinear regime is charac- 
terized by film deformations that are no longer merely 
a linear superposition of contributions with independent 
wave number. Instead, the growth of individual peaks in- 
fiuences the growth of neighboring peaks as determined 
from Eq. (19). The evolving pillars can, for example. 



reposition themselves along directions that are energeti- 
cally favorable in order to maximize the heat flux through 
the air/liquid bilayer and in so doing, can influence the 
in-plane symmetry. This regime can be investigated by 
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holding all remaining parameters fixed while increasing 
Do so as to allow the fluid elongations more time to grow 
before contacting the cooler substrate. This is easily 
achieved in the simulations by either increasing the ac- 
tual plate separation distance, do, or reducing the initial 
film thickness, ho- 

Shown in Fig. ^a) and (b) are two representations of 
the film height H{X,Y,T = Ttop) for AT = 46 °C and 
Do = (a) 3.45 and (b) 7.125. The inset figures depict the 
corresponding FFTs, where the Fourier coefficients have 
been normalized to their peak value and squared for fil- 
tering purposes. The arrow shown has unit length and 
represents the value K'^^^. Contact with the cooler plate 
is achieved at Ttop ~ 1-30 and 1.84, respectively. The 
Fourier transform of the pattern (inset) for the smaller 
value of Do suggests quasi-hexagonal symmetry, with 
some pronounced harmonics in the vicinity of the domi- 
nant peaks. By contrast, the pattern for the larger value 
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FIG. 8: Numerical solutions (top view) of the film thickness, 
H(X,Y,Ttop), from Eq. (191 for different gap ratios. Do = 
do/ho for Ar = 46°C: (a) ho = 100 nm. Do = 3.45, nop = 1.30 
and (b) ho = 40 nm. Do = 7.125, Ttop = 1.84. Inset images 
represent 2D Fourier transforms of corresponding film heights 
viewed from above. Fourier coefficients were normalized to 
the maximum value for each image and squared for improved 
filtering. Unit arrows represent magnitude of most unstable 
wave number, K^-u^^/V^ = 2-k, deri ved from linear stability 
theory [see discussion following Eq. ( 32 1] . 



of Do clearly shows well developed hexagonal symmetry. 
These patterns indicate that the formation of hexagonal 
symmetry is correlated with film depletion near the base 
of nanopillars. For some parameter sets investigated, 
there is also evidence of a bifurcation cascade, in which 
the region halfway in between two adjacent nanopillars 
generates a parasitic protrusion smaller in amplitude but 
similar in shape to the primary nanopillars. This cas- 
cade behavior resembles the dynamics reported in other 
thin film instabilities [S^HSi] , For this cascade to occur, 
the value of Do must be sufficiently large such that the 
growth of the dominant nanopillars consumes a substan- 
tial portion of the interstitial fluid mass. 

In reviewing images of nanopillar formation in the lit- 
erature, it is evident that hexagonal symmetry can oc- 
cur with even small values of Do, as shown in Fig. fljd) 
for which Do = 1.63. If the pillars are allowed to grow 



well beyond the time required for initial contact with the 
cooler substrate, then the dynamics of growth by ther- 
mocapillary stresses will likely continue to draw liquid 
upwards, thereby thickening the diameter of nanopillars 
which bridge the gap in between the two substrates. This 
process will continue to remove film material from the in- 
terstitial regions thereby generating conditions favorable 
to the formation of hexagonal symmetry. In such cases, 
the hexagonal symmetry is likely established well after 
the fastest growing peaks make contact with the cooler 
plate. The mechanism leading to this scenario, however, 
is not included in the model leading to Eq. (14 1. 



2. Effect of substrate tilt 



As discussed in Section II. B, the evolution equation for 



the film height is modified according to Eq. ( 38 1 when the 
confining substrates are subject to a relative tilt. Shown 
in Fig. (|9| are the corresponding results for solutions of 
i?(X, y, Ttop) for the case ho = 100 nm , do — 285 nm 
and AT = 46 °C subject to increasing inclination angle. 
The image shown in Fig. (|9])(c) corresponds to the incli- 
nation angle used in the experiments of Schaffer et al. 




FIG. 9: Numerical solutions (top view) of the film thickness, 
H{X,Y,Ttop), from Eq. (38 1 for different inclination angles 
tan((^) of the cooler substrate: (a) tan(</5) = 4.8 • 10~^ and 
Ttop = 1.23, (b) tan(^) = 4.8 • 10^'' and nop = 1.12, (c) 
tan((^) = 4.8 • 10"^ and nop = 0.86, (d) tan(^) = 4.8 • 10"^ 
and Ttop = 0.47. In all cases, ho — 100 nm, do = 285 nm and 
AT = 46 °C. Inclination angle was imposed along the diag- 
onal of the computational domain such that tan(^)/\/2 =| 
dD/dX 1 = 1 dD/dV \. Schematic diagram indicates that up- 
per right corner (lower left corner) is subject to the smallest 
(largest) wafer separation distance. 

As described in Section II. B, the tilt of the upper sub- 
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strate was defined by the unit vector T|| = (1, l)/\/2- 
The top right corner in the images shown corresponds to 
the region of the film with the smallest gap separation 
distance; likewise, the bottom left corner represents the 
region with the largest gap distance. As such, a lateral 
thermal gradient is established which draws fluid from 
the bottom left corner into the upper right corner. To 
conserve mass in the simulations, fluid exiting the top 
(bottom) boundary was simultaneously replaced by fluid 
entering the right (left) boundary. 



In all the experiments of Schaffer et al. , the confin- 
ing substrates were subject to a relative tilt tan((/3) w 
1 /im/cni. In rescaled units, ta,n{jp) — tan((^)/e where 
e = 'lo/Aj^ax- f'or the experiments in which ho = 100 
nm, do = 285 nm and AT = 46 °C, Ajc^ = 4.8 ^ln\ [Eq. 
(27l], such that tan(^) — 4.8 • 10~^. The corresponding 



tilt angle along the X and Y axes for such experiments 
corresponds to a value of 4.8-10^'^/a/2 w 3.3-10"'^, which 
should lead to the formations observed in Fig. (l9])(c) if 
there were no other considerations or artifacts. 



As evident from the images (b) - (d) , the symmetry of 
the evolving instability transitions from hexagonal-like 
to square-like symmetry due to the lateral bias in ther- 
mal gradient established by the tilt of the cooler sub- 
strate. Even for very small tilt angles, fluid is preferen- 
tially transported toward the upper right corner where it 
accumulates in the form of ridges along the top and right 
boundaries. This accumulation process establishes sec- 
ondary and tertiary parallel ridges spaced apart roughly 
by a distance A^^^^. At longer times, these ridges are 
observed to undergo breakup with a similar lateral spac- 
ing. A relative tilt of the substrates therefore introduces a 
strong lateral bias in thermal gradient which triggers pat- 
tern formation along the domain boundaries instead of 
within the interior, where the instability is generally more 
homogeneously distributed. This specific square symme- 
try observed is therefore a direct consequence of using in- 
clined substrates within a square computational domain. 
Modification of the computational domain shape may al- 
ter the symmetry observed; however, the nanopillars will 
still nucleate along cooler regions of the film. Additional 
studies of the Fourier transforms of the images shown 
in Fig. ([9]) (a) - (d) (not shown) confirm that the fastest 
growing wavelength, Aj^^x; remains unaffected by very 
small tilt angles. In this respect, the measurements of 
Ainax niade by Schaffer et al. with a tilted wafer geome- 
try should not have affected comparison to analytic pre- 
dictions from linear stability theory for films confined by 
parallel wafers. Given the strong influence of edge be- 
havior on the formation of emerging patterns, however, 
care should be taken in experiment to ensure that no ar- 
tifacts, anomalies or asymmetries exist along the edges 
of a film undergoing nanopillar formation if a particular 
array symmetry is desirable. 



V. CONCLUSION 

In this work, we provide evidence that the sponta- 
neous formation of periodic pillar arrays in molten poly- 
mer nanofilms confined within closely spaced substrates 
maintained at different temperatures is due to a ther- 
mocapillary instability. If not mass limited, these pillars 
continue to grow until contact with the cooler substrate 
is achieved. So long as the initial film thickness and sub- 
strate separation distance are sufficiently small that grav- 
itational forces are negligible, there is no critical number 
for onset of instability. In contrast with the conventional 
Benard-Marangoni instability, nanofilms arc prone to for- 
mation of elongations no matter how small the transverse 
thermal gradient. Ultra small gradients, however, lead to 
large values of the most unstable wavelength. In practice, 
very large pillar spacings can be difficult to observe or dif- 
ficult to distinguish from defect mediated bumps which 
also undergo growth from thermocapillary flow. The lin- 
ear stability analysis shows that pillar formations are ex- 
pected in any viscous Newtonian-like nanofilm. Since 
the shear rates are characteristically small, it is expected 
that molten materials of many kinds can be modeled as a 
Newtonian fluid. Pillar arrays formed from polymers like 
PS or PMMA are of commercial interest, however, since 
they solidify rapidly in place once the thermal gradient is 
removed due to their lower glass transition temperatures. 

The analytic results obtained, including the energet- 
ics of nanopillar formation as described by the Lyapunov 
functional, confirm that elongations are caused by the 
predominance of thermocapillary stresses, which far out- 
weigh stabilization by capillary stresses during the later 
stages of development. The increase in thermocapillary 
stresses leads to a rapid decrease in the overall free en- 
ergy of the evolving film. Fourier analysis of the emerging 
structures also indicates a preference for hexagonal pack- 
ing although true hexagonal order cannot be achieved 
if the separation distance is too small since the pillars 
have insufficient time to grow and self-organize before 
making contact with the cooler target. Simulations for 
larger values of Do ~ do/ho show well developed and 
long range hexagonal order. The only limitation of the 
current analysis is the restriction to films of constant vis- 
cosity. While this approximation holds well for simple 
fluids, it is known that the viscosity of polymer melts 
like PS and PMMA exhibit a strong dependence on tem- 
perature. It is therefore expected that fiuid elongations 
undergo an increase in viscosity as the cooler substrate is 
approached. We have examined this effect in detail in a 
separate study [SHj and concluded that while this cooling 
effect slows the growth of pillars, it does not affect the 
pillar spacing in any appreciable way. This is expected 
since the expression for the most unstable wavelength 
given by Eq. (27) is independent of the melt viscosity. 



The linear stability analysis of an initially flat viscous 
film of thickness ho subject only to capillary and thermo- 
capillary forces reveals that the normalized gap spacing 
do/ ho and temperature drop AT strongly affect the value 
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of the most unstable wavelength, AJ^^xj fo^" given material 
parameters. The analysis indicates that the pillar spac- 
ing or array pitch scales as (AT)~^/^, which can therefore 
be tuned in experiment. Direct comparison of Aj^^^ to 
experimental measurements of Schaffer et al. reveals ex- 
cellent agreement with the functional dependence on do, 
namely A^^^ = Ci do^ + C2 do~^''^. The discrepancies 
observed are attributed to a number of factors including 
solvent retention effects in unannealed films and mea- 
surements of the array pitch in vitrified films examined 
after the fluid had experience prolonged contact with the 
cooler substrate. A number of factors not included in 
the model can influence the array pitch since the melt is 
no longer growing in air but migrating and reorganizing 
along the underside of the cooled wafer. 



A linear stability analysis and numerical solutions of 
the nonlinear evolution equation were also conducted for 
a tilted cooler substrate. Such a tilt initially estab- 
lishes both a lateral and vertical thermal gradient. In 
the experiments of Schaffer et al. , the tilt angle was less 
than 0.006°. Numerical simulations of the film height for 
even very small tilt angles confirm that while the dom- 
inant wavelength is unaffected, the in-plane symmetry 
of evolving elongations can transition from hexagonal to 
square-like symmetry. This change is caused by ther- 
mocapillary influx of fluid into the region subject to a 
smaller gap width where the effective surface film tem- 
perature is cooler due to closer proximity with the cooler 
tilted substrate. The elongations in this region also grow 
more rapidly since the effective thermal gradient is larger. 
These results highlight the importance of boundary con- 
ditions in establishing the in-plane symmetry of arrays 
formed as a result of thermocapillary instability in a 
tilted geometry. This observation can also be used to 
advantage to generate large area arrays of different sym- 
metry. 



In conclusion, the results presented here strongly sug- 
gest that thermocapillary stresses play a crucial if not 
dominant role in the formation of pillar arrays in molten 
nanofilms subject to a transverse thermal gradient. Ac- 
cording to the linear stability analysis, nanoscale films for 
which the hydrostatic pressure is completely negligible in 
comparison to capillary and thermocapillary forces will 
promote fluid elongations no matter how small the tem- 
perature difference between the top (cooler) and bottom 
(warmer) substrates. Experiments using lower viscosity 
melts, larger thermal graidents, smaller wafer separation 
distances, and smaller initial film thicknesses should pro- 
duce nanostructures with submicron lateral feature sizes. 
We hope that future studies such as these can assist with 
the design and fabrication of functional devices by tak- 
ing advantage of the inherent regularity, smoothness and 
robustness of self-organized patterns arising from a con- 
trollable hydrodynamic instability. 
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VII. APPENDIX 



To begin, Eq. ( 19 ) is re-expressed in terms of the pa- 



rameter X such that 



dH 



= -Vii • H-" 



nMa 



dr " \2DoH{l + xHY 

which is rearranged according to 



(A-1) 



dH 



-Vii -77^ 



KMa / 1 X X 

^aT \h'i + xh'{i + xH)^ 

1 



3Ca 



Vfi7 . 
(A-2) 



The term proportional to Vyi? is further simplified, 
where 



1 X X 

H l + xH {l + xHY 

In' ^ ^ 1 



Viiiy = 



I+XHJ l+xH 



(A-3) 



By introducing the function ^ = 

{nMa)/{2Do){HH/{l -f xH)] + 1/(1 + xH)) + C„, 
the evolution equation can be recast as 



dH 



= -Vii -H^ 



* 



1 



-y^H 



(A-4) 



where Co is a constant of integration. Eq. ( A-4 ) is then 
multiplied by the quantity ^ = * + Vf.H/{'SCa) to give 



~dH 

*a7 



-*V|| -iJ^Vii*. 



(A-5) 



Since '5 = '^{H), one can apply Leibnitz's rule for differ- 
entiation to find 



dl 



d f"^^'^ 



dr dr 



/ ^{S)dS = ^^ (A-6) 

Jh(t=o) or 



where H{t = 0) = 1 i.e. the initial film is flat and 
uniform. Evaluation of the function / then gives 



I=- 



kMq 

~2d: 



H\n 



H 



1 + xH, 



+ln(l+x) 



-Ci(i/-l)(A-7) 



where Ci denotes a second constant of integration. 



Equation (A-5| is then integrated over the square do- 
ain An = AX AY: 



main Ay = AX Ay 



m 

A., \d^ 3^ 



'' OT 



*V|| ■ H-^X/w-^dXdY, 



(A- 



where ^V|| • (if^Vj|^) can be re-expressed as V|| • 
(t-iJ^Vy*) - i7^(V||*)2. The first term on the right 
hand side vanishes for a fixed domain subject to periodic 
boundary conditions; the integral /^ {dI/dT)dXdY can 

be rewritten as d/dr J^ IdXdY. These simplifications 



can be used to recast Eq. ( A-8 ) into 



^ / I(H)dXdY + i / VlH^dXdY = 

I iJ-'^fv II*) dXdY{A-9) 

A final integration by parts subject to periodic boundary 
conditions simplifies the second integral on the left hand 
side such that 



V\iI^—dXdY = 



1 
1 d 

2d7 



— N\\Hf dXdY 
/ (V||i/)^dA:rfXA-10) 

J A« 
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Equation ( A-9 ) then simplifies to the form 



dr 



I-^N,H? 



QCa 



iJ' 



(V||^) 



dXdY ■ 

~\ 2 



dXdY. 



(A-11) 



Inserting Eq. (A-7) into Eq. (A-11 1 and noting that vol- 



ume conservation within the domain An requires that 
X4 (H - l)dXdY = leads to 



d 



nMa 



^ln(Y^)+ln(l + x) 



~N\\Hf\dXdY = (A-12) 
QCa ^ " ^ J 



if-^ Vii* dXdY. 



Multiplying Eq. (A-12) by the quantity —6 Ca produces 



the final expression for the rate of change of ^, namely 



dT 



ZdXdY 



-6Ca / H'' ( V||* 1 dXdY < 0, 

' A\\ 

(A-13) 



where £ is given by Eq. (34). Since Eq. (A-13) is a non- 



negative quantity, the thin film seeks configurations of 
the interface H in time which minimize 5- 
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